title: GIS Assignment 3: Home Range Data
Figure 1. Gray Wolf (Canis lupus) 
data <- read.csv("15008-15008.csv")
head(data)
##     event.id visible timestamp location.long location.lat external.temperature
## 1 9710323244    TRUE   02:33.0     -112.2716     56.75703                    9
## 2 9710323245    TRUE   00:59.0     -112.2728     56.75678                    5
## 3 9710323246    TRUE   00:32.0     -112.3320     56.77558                    5
## 4 9710323247    TRUE   00:32.0     -112.2713     56.87374                    0
## 5 9710323248    TRUE   01:33.0     -112.3469     56.83778                    0
## 6 9710323249    TRUE   00:58.0     -112.2206     56.80442                    0
##   gps.dop height.raw sensor.type individual.taxon.canonical.name
## 1     2.4     485.46         gps                     Canis lupus
## 2     1.6     476.31         gps                     Canis lupus
## 3     1.2     476.67         gps                     Canis lupus
## 4     1.2     483.20         gps                     Canis lupus
## 5     2.0     481.78         gps                     Canis lupus
## 6     1.8     476.82         gps                     Canis lupus
##   tag.local.identifier individual.local.identifier
## 1                15008                       15008
## 2                15008                       15008
## 3                15008                       15008
## 4                15008                       15008
## 5                15008                       15008
## 6                15008                       15008
##                        study.name utm.easting utm.northing utm.zone
## 1 ABoVE: Boutin Alberta Grey Wolf    422252.2      6291062      12N
## 2 ABoVE: Boutin Alberta Grey Wolf    422180.4      6291035      12N
## 3 ABoVE: Boutin Alberta Grey Wolf    418597.5      6293197      12N
## 4 ABoVE: Boutin Alberta Grey Wolf    422515.1      6304051      12N
## 5 ABoVE: Boutin Alberta Grey Wolf    417826.5      6300137      12N
## 6 ABoVE: Boutin Alberta Grey Wolf    425465.3      6296280      12N
plot <- ggplot() + geom_point(data=data, 
                                   aes(utm.easting,utm.northing,
                                       color=individual.local.identifier)) +
                        labs(x="Easting", y="Northing") +
                        guides(color=guide_legend("Identifier"))
ggplotly(plot)
pal <- colorFactor(
  palette = 'Dark2',
  domain = data$individual.local.identifier
)
leaflet(data) %>% 
  addTiles(group = "Base Map") %>%
  addWMSTiles("https://basemap.nationalmap.gov/arcgis/services/USGSTopo/MapServer/WmsServer", layers = "0", group = "Terrain Map") %>% 
  addCircleMarkers(lng=data$location.long, lat=data$location.lat, radius = .5, label = data$individual.local.identifier, group = data$individual.local.identifier, color = ~pal(individual.local.identifier)) %>%
  addLayersControl(
    baseGroups = c("Base Map", "Terrain Map", "ESRI"), overlayGroups = c(unique(data$individual.local.identifier)),
    options = layersControlOptions(collapsed = FALSE)) %>%
    addScaleBar() %>%
    addMiniMap(zoomLevelOffset = -4, position = 'bottomleft')
lapply(split(data, data$individual.local.identifier), 
       function(x)write.csv(x, file = paste(x$individual.local.identifier[1],".csv", sep = ""), row.names = FALSE))
## $`15008`
## NULL
## 
## $`35260`
## NULL
files <- list.files(path = ".", pattern = "[0-9]", full.names = TRUE)
utm_points <- cbind(data$utm.easting, data$utm.northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=12 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.01, min(proj_lat.lon$x)-0.01), 
                  c(min(proj_lat.lon$y)-0.01, max(proj_lat.lon$x)+0.01), 
                  type = "bing")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=12 +datum=WGS84")

autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "purple", fill=NA, size=1)) +
  geom_point(data=data, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 1, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "red", fill=NA, size=1)) +
  geom_point(data=data, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

library(maptools)
if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit() #This resolved the following error message from pbapply: "Error: isTRUE(gpclibPermitStatus()) is not TRUE"
## [1] TRUE
mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual.local.identifier)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "blue", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=long, y=lat), alpha=0.8) +
    geom_point(data=mcp.points, aes(x=x, y=y)) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)
## [[1]]

## 
## [[2]]

## 
## [[3]]

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 95)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual.local.identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.8) +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)
## [[1]]
## Warning: Use of `kde.poly$long` is discouraged. Use `long` instead.
## Warning: Use of `kde.poly$lat` is discouraged. Use `lat` instead.

## 
## [[2]]
## Warning: Use of `kde.poly$long` is discouraged. Use `long` instead.

## Warning: Use of `kde.poly$lat` is discouraged. Use `lat` instead.

## 
## [[3]]
## Warning: Use of `kde.poly$long` is discouraged. Use `long` instead.

## Warning: Use of `kde.poly$lat` is discouraged. Use `lat` instead.

Figure 2. Gray Wolf Pups